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ABSTRACT 


Numerical  results  are  presented  for  the 
diffusive  growth  in  time  of  a  barium  cloud  in 
a  highly  (but  not  infinitely)  conducting  iono¬ 
sphere.  Polarization  effects  produce  steepening 
of  the  backside  and  a  thinning  and  indentation  on 
the  front  side.  A  crude  estimate  of  the  time  re¬ 
quired  for  appreciable  distortion  is  given. 

Numerical  results  are  also  given  for  the 
linearized  equations  of  motion  about  an  ellipsoi¬ 
dal  gaussian  equilibrium  distribution.  A  new 
numerical  method  yields  the  fastest  growing  eigen¬ 
value  and  eigenfunction  and  has  been  pushed  to 
49  x  49  size  grids  covering  the  cloud.  The  resul¬ 
tant  eigenfunctions  (in  the  plane  perpendicular  to 
the  B  field)  are  concentrated  in  the  rear  of  the 
cloud  and  are  elongated  parallel  to  the  cloud  drift 
velocity.  Finally,  it  is  shown  that  results  already 
available  in  the  original  E  X  B  instability  paper 
can  be  used  to  obtain  an  estimate  of  the  striation 
scale  size  and  onset  time. 
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I. 


INTRODUCTION 


We  have  previously  reported^)  the  development  of  an 
analytic  model  for  the  slow  diffusive  growth  in  time  of  a 
barium  cloud.  This  model  assumed  a  background  ionosphere 
which  was  highly  conducting  but  no  longer  infinitely  so. 

This  generalization  of  previous  work^)  allowed  polarization 
charges  to  accumulate  on  the  cloud  and  thus  permitted  the 
cloud's  shape  to  distort  away  from  a  gaussian  ellipsoid. 

In  Section  II  of  this  report,  we  give  details  of  numerical 
calculations  of  this  effect  and  present  some  illustrations 
of  it.  We  are  also  able  to  make  a  crude  estimate  of  the 
time  required  for  appreciable  steepening  to  occur.  This 
is  described  in  Section  III. 

In  an  earlier  technical  report^)  We  also  described 
numerical  solutions  of  the  linearized  equations  of  motion 
of  the  cloud.  These  equations  were  linearized  about  the 
equilibrium  gaussian  ellipsoid  derived  in  ref.  2.  We  also 
noted  in  this  report  that  the  particular  numerical  scheme 
used  was  rather  inefficient  and  had  been  pushed  close  to 
its  practical  limit  by  the  time  the  cloud  had  been  repre¬ 
sented  by  a  9  x  9  grid.  We  have  now  changed  to  an  entirely 
different  and  much  more  efficient  scheme.  Details  and 
illustrations  of  the  results  are  given  in  Sect.  IV.  Finally, 
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Finally,  we  note  that  it  is  possible  to  make  crude  esti¬ 
mate  of  the  scale  size  of  striations  and  estimate  their 
onset  time  by  a  simple  extrapolation  of  results  already 
in  the  original  1963  paper  on  the  E  X  B  instability. ^ 
Details  are  given  in  Sect.  V. 


II.  NUMERICAL  CALCULATION  OF  CLOUD  DISTORTION 


la  the  previous  semi-annual (^),  (hereafter  referred  to 
as  I),  an  equation  was  derived  for  the  growth  by  diffusion 
and  convection  of  a  cloud  immersed  in  a  highly  conducting 
(but  not  infinitely  conducting)  background  ionosphere. 

We  also  noted  some  typical  numerical  results.  The  pur¬ 
pose  of  this  section  is  to  further  describe  the  numerical 
scheme  used. 

The  basic  equation  is  equation  (17)  of  I. 
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and  £.  ,  the  first  order  correction  (in  the  ratio  of  con- 

ductivities)  to  the  electric  field  is  given  by  equation  (14) 
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where 


i  ‘  ~  t.  f  V'4  . 


The  initial  state  is  assumed  to  have  a  Gaussian  profile 
of  width  a,  and  is  normalized;  to  the  total  number  of  particles 
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It-d?-  (3) 

To  proceed  further,  we  must  fix  a  coordinate  system. 

We  assume  that  the  external  electric  field  points  in  the 
positive^/  direction,  the  magnetic  field  in  the  positive  2: 
direction,  and  hence  the  E  X  B  drift  is  in  the  positive  X 
direction.  (We  remind  the  reader  that  the  equations  are 
expressed  in  coordinates  which  move  with  respect  to  the  rest 
frame  of  the  neutrals  with  velocity  equal  to  the' O’th  order 


drift  given,  by 


r  -(V-/  )  &  *  k 


cf.  equation  (12)'  of  I) 

We  also  reduce  the  equation  to  nondimens ional  form  by 
measuring  length  in  units  of  the  initial  Gaussian  radius 
and  time  in  units  pf  : 

Jl  -  -22  / «_ 

r-  r  (4) 

Substituting  (2)  into  (1),  specializing  to  our  coordinate 
system  and  using  (4),  we  obtain  our  basic  equation 
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In  the  above  *we  have  neglected  terms  of  order  / Q^t) 


so  that 


and 


c 

B 


In  addition,  we  have  defined  the  two  parameters; 

T  -  g  ~S>° 

which  is  the  ratio  of  the  cross-cloud  conductivity  to  that 
of  the  background  and 


y  —  g  ^ 
d  "  *Cr<iK) 

which  is  the  normalized  external  electric  field  strength. 
The  solution  therefore  depends  only  upon  the  specification 
of  the  three  parameters  (/hT).  >  A  and  Y . 

Numerical  solutions  to  equation  (5)  have  been  obtained 


using  the  alternating  direction  implicit  scheme  (A.D.I). 

In  this  method  the  differential  equation  (5)  is  replaced  by 
the  following  two  coupled  finite  difference  equations. 


where  A)  refers  to  the  time  step  and  c  and ^  to  the  grid 
point  in  X  and Jj  respectively. 

Notice  that  in  the  first  equation  (7a)  the  differences 
in  H  are  advanced  one  time  step,  whereas  in  the  second  the 
differences  in  7  are  advanced  one  step.  (The  coefficients 
A,  B,  and  C  are  evaluated  "in  the  middle",  that  is  they 
are  evaluated  at  t.  **  (n+1)  for  both  equations.) 

This  method  can  be  shown  to  be  much  more  stable  than 
explicit  schemes  in  which  all  spatial  differences  are  taken 
at  the  present,  known,  time  step.  (For  the  case  X  “  0, 
A.D.I,  gives  an  unconditionally  stable  solution).  On  the 
other  hand,  the  method  is  computationally  much  simpler  than 
a  completely  implicit  scheme  in  which  the  differences  in 
both'  X  and  are  advanced  at  the  same  time.  The  reason, 
as  can  readily  be  seen  from  equation  (7),  is  that  the  set 
of  equations  to  be  solved  at  each  time  step  are  tridiagonal 
in  form  and  are  thus  susceptible  to  solution  by  a  very 
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economical  algorithm.  We  have  used  a  double  precision 
subroutine  which  employs  Gaussian  elimination  with  partial 
pivoting. 

We  set  up  the  equations  on  a  square  grid,  choosing 

.20834,  which  proved  sufficiently  small  to  guar¬ 
antee  accuracy,  and  A T  as  large  as  possible  while  still 
ensuring  stability.  Some  error  was  inevitably  introduced 
by  the  use  of  a  finite  sized  grid  for  a  problem  which  is 
defined  on  an  infinite  plane.  We  adopted  the  scheme  of 
setting,  all  of  the  grid  boundary  points  to  zero  at  each  time 
step.  For  several  of  the  runs  in  which  the  cloud  diffused 
rapidly i  we  used  an  expanding  grid  scheme,  that  is  as  time 
developed  and  the  cloud  grew,  we  increased  the  overall  size 
of  the  grid.  In  all  cases  the  value  of  the  density  at  the 
boundary  points  was  no  greater  than  five  orders  of  magnitude 
less  than  the  maximum  value  of  the  density. 

In  general,  errors  arose  from  three  sources: 

1. )  truncation  error  due  to  the  finite  grid  spacing, 

2. )  machine  round-off  error,  and 

3. )  error  due  to  the  finite  size  of  the  grid. 

The  first  was  checked  by  choosing  different  values  of  grid 
spacing;  the  second  by  using  two  compilers,  one  which  rounds 
off  to  the  next  larger  integer  and  the  other  to  the  smaller; 
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and  the  third  by  varying  the  grid  size;  The  resulting  errors 
from  all  causes  are  less  than  10%. 

Some  computer  results  are  shown  in  Figs.  1-5.  In  all 
cases,  we  set  A  *  0.1,  $'m  20  and  varied  QT.7~j)^  .  The  grid 

scale  in  these  plots  is  in  kilometers  with  4.8  km.  To 
get  some  idea  of  the  real  time  corresponding  to  the  various 
values  of  fused,  note  that  Jh'z  .045  km/&-  (Secede  1,  Dogwood) 
and  thus  t  ('real]  »  (f- r~  127Tseconds.  Note  the  thin¬ 
ning  and  indentation  on  the  front  side  of  the  cloud  and  the 
steepening  in  the  rear.  The  arrow  marked  ^represents  the 
actual  direction  of  drift  of  the  cloud  relative  to  the  neutrals. 
Note  the  symmetrical  distortion  of  the  cloud  in  the  large 
limit.  This  symmetry  is  readily  verified  by  examination  of 
the  coefficients  in  Eq.  (6)  in  this  limit.  The  leading  terms 
in  B  are  proportional  to  /and  while  those  in  C  are  pro¬ 


portional  to  xy  . 


III.  ESTIMATE  OF  STEEPENING  TIME 

One  may  also  obtain  a  crude  estimate  of  the  "steepening 
time"  by  examining  the  coefficients  in  Eq..  (6).  In  the  large 
(h r),  limit,  the  leading  terms  are  in  coefficients  B  and  C 


and  are  of  order 


or  ir  A  ft  (flry 


Hence  the  distortion  time  7^  is 


k  A  ¥  ^17*3+- 

This  is  in  units  of  the  diffusion  time  In  rea* 

time*  we  have 
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IV.  NEW  NUMERICAL  RESULTS  FOR  THE  LINEAR  INSTABILITY 
EIGENFUNCTIONS 

A:  Linear  Equation: 

Some  preliminary  numerical  calculations  of  the 
fastest  growing  linear  eigenvalues  were  presented 
in  an  earlier  semiannual  report.  ^  We  now  have 
changed  to  a  much  improved  numerical  scheme  and 
report  considerably  changed  results  on  much  larger 
grids.  The  starting  point  for  this  calculation  is 
the  same  set  of  equations  as  shown  in  Eq.  (2.10) 
of  ref.  3  above.  These  are 

[w  t  Jj-X-  K  I  (a/  x') 

7/1  Z  ■  (  A  x  tX  m  +  1  <e)  l ±x  X 

xfc  A?  -  Ac  •  Pj.  N-  L  £t,x&  ■  Z.A/J 
-  (ar\  pi  i'  [a  Fjl  '  ?a7  -t  p  a7  X  n  <(■] 

+  L~  ft  *-/<-]  €•  V(*X) 

+  f~p-i)  +  Cv-tX  pi]  f,  6  (»*)  —  o 


(10) 
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where  we  have 

replaced 

— >  -  /S 

and  redefined 

Sj. 

K  H 

%/y 

L  s 

Otherwise,  all  symbols  have  the  same  significance  as  in 
ref.  3.  Next  we  reduce  this  equation  to  dimensionless  form 
by  defining 

£  = 

Ckj  ~sL  La  XT 


s  /A 

*■  2  - 

*(  Ti  +  r.) 

y—  ~  V 

'  >- 

and  further  simplify  by.  dropping  terms  of  order  (EIt)+  //XlT'j 
(consistent  with  our  expansion  in  powers  of  /  z);/  \  . 

The  resultant  equation  in  Cartesian  coordinates,  with  EL  0 
in  the  positive  y  -direction,  is: 

D-x  3JL 

^  jf  4  F~Tj  +**£<<  = 

-y-  an 


where 
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We  take  /\/  to  have  a  gaussian  radius  ^  at  the  origin  of  time. 


Thus 


B:  Numerical  Scheme 


^  Kfi  /-  'J/  (/  +  r)J 


The  two  coupled  equations  above  can  be  written  in  the 


form: 


l;  y.  t  t-t  <*  -  ^  x 


(13a) 


(13b) 
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If  one  replaces  the  differential  operators  by  the  finite 
difference  approximation, 

^  UL  ilu 

iv  Z-v'  (i4) 

;>  aX  '  y 

ThT  ~  £oc rj 

cr  a  4  k 

f  £ 

the  operators  Z  and  become  matrix  operators.  The  set 
of  equations  are  closed  by  specifying  homogeneous  boundary 
conditions  on  X and  P  •  If  we  now  subtract  Eq.  (13b)  from 
(13a)  we  obtain 

(Zr-c)*  -  -<v-o*  (15) 

Assuming  is  nonsingular,  we  can  invert  the  matrix  and 

substitute  the  resulting  equation  for  (p  back-  into  (13a). 

We  can  then  solve  (13a)  for  the  eigenvalue  and  eigenfunction 


by  any  of  a  number  of  algorithms. 

This  straightforward  procedure  was  the  method  originally 
employed  for  this  problem^).  Unfortunately,  it  suffers 
from  the  limitation  of  being  very  uneconomical.  To  see  this 
we  observe  that  an  grid  has  /^points.  (When  we  refer  to 
the  dimensions  of  a  grid,  we  do  not  consider  the  boundary 
points  which  do  not  enter  the  matrix  equations  due  to  the 
homogenous  boundary  conditions  onXand<?).  Eq.  (13)  includes 
one  equation  for  each  point  and  hence  the  coefficient  matrices 
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will  be  arrays  of  dimension />t  by  fn  for  a  total  of  /w  elements. 
The  storage  requirements  of  the  computer  therefore  increases 
as  the  fourth  power  of  the  grid  dimension  and  the  method 
very  quickly  outgrows  the  capacity  of  the  computer.  Indeed 
it  was  found to  be  impossible  to  consider  grids  with  more 
than  11  x  11  points.  This  is  hardly  a  sufficiently  fine  mesh 
to  discern  the  complicated  structure  of  an  unstable  mode. 

To  obtain  a  more  economical  procedure,  we  first  observe 

that  most  of  the  elements  in  the  matrices  of  coefficients 

are  zero.  Indeed,  it  is  clear  from  Eq.  (14)  that  there  are 

only  5  possible  ijon-zero  bands  of  elements  running  diagonally 

2. 

across  each  coefficient  matrix;  a  total  of  at  most  only -57*1  non¬ 
zero  elements.  The  problem  enters  in  the  solution  of  Eq.  (15) 
since  the  inversion  of  the  matrix  (4^~ZX)  will  in  general  des¬ 
troy  the  band  structure  and  produce  a  full  matrix. 

Our  solution  to  this  problem  involved  a  double  iterative 
procedure.*  We  solved  equation  (13a)  for  the  eigenvalue  and 
eigenvector  ^  iteratively  via  the  power  method^),  obtaining 
for  each  iterate  by  solving  the  set  of  equations  (15).  That 
is,  an  initial  approximation  for  ^  was  randomly  generated  and 
substituted  into  the  R.H.S.  of  equation  (15).  The  resulting 
set  of  inhomogeneous  linear  equations  was  solved  for  Cp  using 
point  successive  overrelation^) ,  a  variation  of  the  usual 
Gauss-Seidel  with  a  parameter  to  accelerate  convergence. 
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The  result  for  p  was  substituted  into  (l3aj  and  a  new  iterate 
for  %  was  obtained  using  the  power  method.  This  was  substi¬ 
tuted  back  into  the  R.H.S.  of  Eq.  (15)  and  then  the  process 
was  continued. 

This  process  in  general  converges  to  the  eigenvector 
corresponding  to  the  eigenvalue  with  largest  absolute  magni¬ 
tude^).  However/,  we  desire  the  fastest  growing  eigenmode', 
that  £s  the  eigenmode  whose  eigenvalue  has  the  largest  real 
part.  To  obtain  the  latter  eigenmode,  we  solved,  in  place 
of  the  original  equation  which  takes  the  form 

A  'j'  ~  40 

(4  is  a  matrix  operator),  the  equation': 

(A  +  px)X  =  *  X  (16) 

where  f>  is  a  real  constant  and  X  is  the  unit  matrix.  It  is 
clear  that  if  we  kept  all  of  the  original  eigenvalues  in  the 
complex  plane,  the  effect  of  adding  the  fL  term  in  Eq.  (16) 
will  be  to  shift  all  of  the  values  parallel  to  the  real  axis 
by  an  amount  -f*  .  One  can  readily  show  by  geometrical  argu¬ 
ments,  that  there  always  exists  a  value  for  the  shift 
which  shifts  the  spectrum  of  eigenvalue  such  that  the  original 
eigenvalue  (or  complex  conjugate  pair  of  eigenvalues)  with  the 
largest  real  part  becomes,  after  the  shift,  the  one  with  the 
largest  absolute  value.  The  process  will  then  converge  to 
this  eigenmode,  and  one  can  substract  the  shifty  to  obtain 
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the  original  value  for  the  eigenvalue. 

C:  Results : 

We  have  run  a  series  of  numerical  solutions  for  grid 
sizes  13  x  13,  19  x  19,  25  x  25,  37  x  37  and  49  x  49.  For 
the  purpose  of  these  runs  we  assumed  the  following  .parameter 
values : 

y  ~  3  " 

-  So 

h-  - 

n. 

We  also  assumed  that  the  grid  covered  three  ^.-folding 
lengths  of  the  equilibrium  density,  i.e.  we  assumed  ~ 

where  /?0  is  the  distance  from  the  center  of  the  grid  along  pc 
^y)  to  the  edge. 

Figures  6-9  show  plots  of  the  density  (07  =•//£)  for  the 
four  grids:  13  x  13,  19  x  19,  25  x  25,  and  37  x  37  for  the 
case  T  -  (the  49  x  49  grid  converged  too  slowly  for  a 

reasonable  result  for  the  eigenvector;  we  have  used  only  the 
result  for  the  eigenvalue  -  see  below).  Notice  that  all  four 
grids  show  striations  along  the  lower  left-hand  corner  of  the 
cloud.  The  amplitude  of  the  perturbation  density  throughout 
the  rest  of  the  cloud  is  smaller  by  at  least  three  orders  of 
magnitude  compared  with  the  maximum  value  in  the  striations. 
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Moreover,  the  striations  are  directed  approximately  parallel 
to  the  equilibrium  drift  velocity:  S0i 

Hence,  they  point  out  the  rear  end  of  the  cloud,  a  result  in 
accord  with  observation. 

In  Table  1  we  list  the  eigenvalues  corresponding  to  the 
eigenfunctions  shown  in  figure  6-9';  We  observe  that  the 
eigenvalues  are  complex;  but  their  real  parts  are  of  order 
ten  times  their  imaginary  .parts. 

The  question  of  convergence  to  a  unique  eigenvalue  is 
a  matter  of  continuing  concern.  In  figure  10  we  plot  the 
real  part  of  the  eigenvalue  as  a  function  tof  grid  spacing  AX. 
One  would  hope  that  as  av  —»>  c  ,  the  value  of  A  would  level 
off  and  achieve  a  constant  value.  Unfortunately,  as  can  be 
seen  from  figure  10,  even  for  a  49  x  49  grid  this  has  not 
occurred.  We  therefore  conclude  that  as  we  go  to  finer  grids 
we  are  picking  up  even  faster  growing  modes.  To  test  whether 
this  process  converges,  we  plotted  log  4/  as  a  function  of  A 
If  the  curve  were  concave  downward  the  process  would  converge 
whereas  if  it  were  concave  upwards  it  would  not.  The  results, 
displayed  in  figure  11  show  the  curve  to  be  very  slightly 
concave  downward  although  the  result  does  strain  somewhat 
the  allowed  errors  placed  upon  the  data  points. 


V.  Approximate  Model  for  Estimating  Striation  Scale 

Size  and  Onset  Time. 

If  one  is  willing  to  settle  for  use  of  a  slab  model 
to  represent  the  cloud,  it  should  be  possible  to  use  the 
results  already  present  in  Ref.  4.  We  shall  show  how  to 
do  this  using  simple  approximations.  All  notation  will 
follow  that  paper  (including^  being  the  magnetic  field 
direction! ) 

The  basic  result  in  Ref.  4  is  contained  in  Eq.  (3.18). 

Let  us  now  simplify  by  assuming  that  the  most  unstable  mode 
is  one  with  very  long  wavelength  in  the  magnetic  field  direc¬ 
tion.  This  is  because^ -variation  leads  only  to  damping. 

Hence  letting  A  c  our  result  reduces  to: 

r  at  /\  p 2 )  -  ®  A  )  J 

(a)  —  ..  - _ — *  — ■  ■  * 


where 


r  -  -  Ay  (/i  Vu  A.  ) 

®  =  -  e0  mt  [M,  •'  *  r4J  f 


Ax=  &/*)*. 


where  we  already  assumed  >>  / .  It  is  obvious  that  the 
second  term  in  the  denominator  is  small  compared  to  the  first 
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for  >  I  .  Note  that 


fijsv la  \ 

("TC77 


(rrMt 

L  X 


A 


and 


^  *  /l' 


Hence  we  drop  this  second  term  and  also  neglect  /Uj_  compared 

to  /<J~, 

Let  us  compare  the  three  terms  in  <&  .  The  ratio  of  the 


first  term  to  the  second  is: 

f- 


O'  Fo  Mc  /0 


oi  )  fr/"'- 

We  assume  that  this  is  large  compared  to  unity  (in  most  cases, 
it  is  of  order  20  -  100).  Hence,  we  neglect  the  second  term. 
By  the  same  means,  the  ratio  of  the  first  to  the  third  is 


—  —  .i. 

^  tc  Ma  A 


>.  ,  i. 

>  A  Y  >  / 


Hence 


&  *  -  F„  />*>«/  «/3 

The  resulting  expression  for  the  growth  rate  is: 


CO.  . 


J-  ‘ilk. 


%(at)+  /\r  v0  .  (i7) 

The  first  term  is  recognized  as  the  damping  due  to  ambipolar 
diffusion  in  the  perpendicular  direction  (the  result  of  taking 
the  second  term  is  destabilizing  if  tc  is  positive 
(the  usual  necessary  condition  for  the  E  x  B  instability). 
There  is  a  critical  wavelength  for  given  i~0.  To  see  this, 
note  that  the  damping  term  is  finite  for  0  and  goes  to 
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infinity  as  'v’  ,  The  destabilizing  term  goes  to  zero 

as  and  goes  to  a  constant  value  as  .  A 

conceptual  plot  of  these  terms  (for  small  /i,  )  is  shown  below. 


the  destabilizing  .term  increases  until  the  curves  cross. 

/  z.  _ 

This  gives  the  threshold  values  of  and  L'o  .  From  the 
shape  of  the  curves  and  the  fact  that  both  have  their  turnover 
in  the  vicinity  of  ^  it  is  obvious  that  the  threshold 

value  of  will  be  in  the  vicinity  of  . 

One  can  also  view  the  onset  of  instability  as  a  function 
of  the  growth  in  cloud  size  while  holding  f0  fixed.  Growth 
means  that  £  increases  just  as  fa  does.  Thus  as  the 

cloud  grows,  the  stabilizing  term  decreases  at  its  lower  end, 
while  the  destabilizing  term  increases  at  the  lower  end  and 
decreases  at  the  upper  end.  Again  when  the  two  terms  cross, 
we  have  the  critical  wavelength  and  the  critical  size  of  the 
cloud  (which  can  be  related  to  an  onset  time  for  any  given 
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theory  or  empirical  set  of  values  of  cloud  size  as  a  function 
of  time). 

Note  that  in  the  derivation  in  Ref.  4,  the  length^  repre¬ 
sents  a  distance  in  the  direction  of  the  density  gradient  over 
which  the  perturbation  extends  in  that  direction.  It  is  assumed 
small  compared  to  the  density  gradient  itself.  It  may  be  possi¬ 
ble  to  get  useful  estimates  of  scale  size  and  onset  times  by 


arbitrarily  assuming 


f)  ~  _L 
st  -  r~  / 


clMo  \~  t 


<r  ^ 


For  example,  if  we  set  ~(^)and  use  the  estimate  above,  it 
is  easy  to  calculate  the  critical  scale  size  from  Eq.  (17). 
Equaling  the  two  terms  on  the  r.h.s.,  one  has 


where 


Mf)  (ti-  M  -  c-^- 

V*'  r-J  3 


M  - 


_  /_L_  A  ! 


Hence 


where 


<*-o  ^  K) 

^  /££  771  A-  , 

(lim  tt) 
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If  we  take  typical  values: 


we  obtain 


Pjp  Jif  /»%«<  . 

^  '>  So  /wn  / 

n/r.  ^  'A- 


/  s~  x/e 


A  very  crude  estimate,  of  course,  but  not  entirely  unreasonable 
since  C/)r]_  ~  /C?  . 
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FIG.  I  DENSITY  CONTOURS  FOR  (QT)+  =50  AT  T  =  .I875 
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FIG.  2  DENSITY  CONTOURS  FOR  (&T)+  =  50  AT  T  =  . 5615 


FIG.  3  DENSITY  CONTOURS  FOR  (&T)+  =  !Q.O  AT  T=  1.875 
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FIG.  4  DENSITY  CONTOURS  FOR  (&T)+=I0.0  AT  T  =4.6875 


FIG.  5  DENSITY  CONTOURS  FOR  (&T)+=2.0  AT  T=14.06 
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FIG.  6  DENSITY  PERTURBATIONS  13x13  GRID 
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DENSITY  PERTURBATIONS  19x19  GRID 


DENSITY  PERTURBATIONS  25x25  GRID 


THE  REAL  PART  OF  X  AS  A  FUNCTION 
GRID  SPACING  AX 


^-49x49  GRID 


13x13  GRID 


LOG  (AX)  AS  A  FUNCTION  OF  THE  REAL  PART 
OF  X 


TABLE  1 


List  of  Eigenvalues 


&Cul) 


versus 


Grid  Size 
1 3x  1 3 
19x19 
25x25 
57x37 


Grid  Spacinp 
1/2 
1/5 
3/4 
1/6 
3/8 


Grid  Size 

Eigenvalue 
4  7 . 5  ±  6 . 4  i 
61.2±3.2i 
7 1 . 2  ±  4 . 9  i 
78. 5  ±  8  i 


49x4  9 


85 . 8±1 2i 
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